In [1]:
pwd
Out[1]:
In [2]:
cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
In [3]:
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar
import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [4]:
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
In [5]:
from scipy.stats import boxcox, probplot
In [6]:
bios = gpd.read_file('biomes_continents')
In [7]:
bios.columns
Out[7]:
In [9]:
area = bios['SUM_Area']*1e-6
In [12]:
data_columns =bios.columns[4:4+17*12]
In [13]:
inter_cols = {}
for yr in range(17):
inter_cols[yr] = data_columns[12*yr:12*yr +12]
In [14]:
intra_cols = {}
for mo in range(12):
intra_cols[mo] = [data_columns[mo + yr *12 ] for yr in range(17)]
In [15]:
for yr in range(17):
inter = bios[ inter_cols[yr]].sum(axis=1)/area
label = 'inter_' + str(yr +2001)
bios= bios.assign(label=inter)
bios = bios.rename(columns={'label': label})
In [16]:
inter_columns = bios.columns[-17:]
In [17]:
inter_max = bios[inter_columns].max().max()
In [18]:
for mo in range(12):
intra = bios[ intra_cols[mo]].sum(axis=1)/(17*area)
label = 'intra_{:02}'.format(mo)
bios = bios.assign(label=intra)
bios = bios.rename(columns={'label': label})
In [19]:
intra_columns = bios.columns[-12:]
In [20]:
intra_max = bios[intra_columns].max().max()
In [21]:
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100):
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_aspect('equal')
ax = df.plot(ax =ax, column=column, cmap='OrRd',vmin=vmin, vmax=vmax)
ax.set_axis_off()
norm = cls.Normalize(vmin, vmax)
cmmapable = cm.ScalarMappable(norm, 'OrRd')
cmmapable._A = []
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig(name, dpi=dpi, bbox_inches='tight')
Inter-anual¶
In [22]:
lambdas = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column, lbd = boxcox(bios[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
In [23]:
lambdas
Out[23]:
In [24]:
probplot(bios.inter_2017, plot=plt.gca());
In [25]:
probplot(test_column, plot=plt.gca());
In [26]:
maxs = []
mins = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
In [27]:
box_min = min(mins)
box_delta =max(maxs) - box_min
In [28]:
for yr in range(17):
column = 'inter_{}'.format(yr +2001)
test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
bios = bios.assign(test_column =test_column)
yr_name = str(yr +2001)
#"Burnt area pixels (#/km2) by biome in 2001
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in ' +yr_name
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'inter_biomes_parts_index_{}'.format(yr +2001)
draw_map(bios,'test_column',title=title,name=name)
Intra-anual¶
In [29]:
lambdas = []
for mo in range(12):
column= 'intra_{:02}'.format(mo)
test_column, lbd = boxcox(bios[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
In [30]:
lambdas
Out[30]:
In [31]:
probplot(bios.intra_11, plot=plt.gca());
In [32]:
probplot(test_column, plot=plt.gca());
In [33]:
maxs = []
mins = []
for mo in range(12):
column= column= 'intra_{:02}'.format(mo)
test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
In [34]:
for mo in range(12):
column = 'intra_{:02}'.format(mo)
test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
bios= bios.assign(test_column =test_column)
mo_name = calendar.month_name[mo + 1]
#"Burnt area pixels (#/km2) by biome in Janury (2001-2017)"
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in ' +mo_name + ' (2001-2017)'
#title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by climate 2001-2017 ' + mo_name
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'intra_biomes_parts_index_{:02}'.format(mo)
draw_map(bios,'test_column',title=title,name=name)
Total¶
In [35]:
total = bios[data_columns].sum(axis=1)/17
total = total/area
bios= bios.assign(total=total)
In [38]:
column = 'total'
test_column, lbd = boxcox(bios[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
bios = bios.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in (2001-2017)'
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'total_biomes_parts_index'
draw_map(bios,'test_column',title=title,name=name)
